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Abstract 

The cosmic ray source spectrum produced by AGN (Active Galactic Nu- 
cleus) jets is calculated. A distinctive feature of these calculations is the 
account for the jet distribution on kinetic energy. The expected cosmic ray 
spectrum at the Earth is determined with the use of a simple numerical code 
which takes into account interactions of ultra-high energy protons and nuclei 
with the background radiation in an expanding universe. 

Keywords: 

ultra-high-energy cosmic rays; acceleration; propagation; AGN jets 
1. Introduction 

The origin of cosmic rays with energies E > 10^^ eV remains a key prob- 
lem of cosmic ray astrophysics. For a uniform distribution of sources at the 
Hubble scale c/i^o ~ 3 x 10^h~^ Mpc, a characteristic cutoff in cosmic ray 
spectrum at 5 x 10^^ eV is developed (here c is the velocity of light, Ho = lOO/i 
km(s Mpc)~^ is the Hubble parameter at the present epoch, h = 0.7). The 
cutoff arises because energetic protons lose energy by electron-positron pair 
and pion production (the GZK effect by Greisen (1966) and Zatsepin & 
Kuz'min (1966)) and energetic nuclei in addition are subject to photodisin- 
tegration (Puget et al., 1976). The observed suppression of the cosmic ray 
flux at energies above ~ 5 x 10^^ eV (Abbasi et al., 2008; Abraham et al.. 
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2010) confirms the presence of tlie GZK cutoff although the suppression due 
to acceleration limits in cosmic ray sources can not be excluded. A debated 
topic is the value of energy Ec in the range 10^"^ eV to 10^^ eV where the 
Galactic component gives way to the extragalactic component in the energy 
spectrum of cosmic rays observed at the Earth, see AUard (2009); Berezinsky 
(2009) for discussion. 

The present knowledge about the highest energy cosmic rays was mainly 
acquired from the two major experiments HiRes (Abbasi et al., 2009; Sokol- 
sky et al., 2010) and Auger (Abreu et al., 2010; Abraham et al, 2010). 
The mass composition of the highest energy cosmic rays remains uncertain. 
The interpretation of HiRes data favor the proton composition at energies 
10^^ — 5 X 10^^ eV, whereas the Auger data indicate that the cosmic ray 
composition is becoming heavier with energies changing from predominantly 
proton at 10^*^ eV to more heavy composition and probably reaching the 
pure Iron composition at about 5 x 10^^ eV. The mass composition interpre- 
tation of the measured quantities depends on the assumed hadronic model 
which is based on not well determined extrapolation of the physics from 
lower energies. The found angular distributions of events exceeding 6 x 10^^ 
eV (these particles can reach the Earth from the distances less than ~ 100 
Mpc) were also different in these two experiments: the HiRes data showed 
the isotropic distribution; the Auger data demonstrated the correlation with 
AGN or objects having a spatial distribution similar to the distribution of 
matter in the nearby Universe with the largest excess from the region of the 
sky around the radiogalaxy CenA. The first data on the spectrum, compo- 
sition and anisotropy from the new Telescope Array experiment support the 
HiRes results (Thomson, 2010). 

The list of potential sources which could give the observed flux of the 
highest energy cosmic rays includes AGNs, gamma-ray bursts, magnetars, 
interacting galaxies, large-scale structure formation shocks and other objects, 
see e.g. review by Torres & Anchordoqui (2004). 

It is assumed below that the AGN jets are the main extragalactic sources 
of ultra-high energy cosmic rays. The kinetic energy of an individual jet de- 
termines the maximum particle energy and the cosmic-ray power it produces. 
The cases of power-law and delta-shaped spectra of cosmic rays produced by 
an individual jet are considered. The distribution of AGN jets over the ki- 
netic energy shapes the average source spectrum of accelerated particles. A 
simple numerical code is used to calculate the expected intensity of cosmic 
rays at the Earth. 
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2. Transport equations for cosmic rays 

We start with a consideration of cosmic ray propagation in the intergalac- 
tic space. The approach based on the equation for cosmic ray density is 
appropriate for the purpose of the present research. 

Let us assume a uniform distribution of cosmic ray sources in the universe. 
The cosmic ray number density per unit energy N{E, t) and the integral num- 
ber density N{> E,t) = dEN{E,t) are then also uniformly distributed, 
i.e. does not depend on position r. The total number of energetic particles 
in the comoving volume V is conserved and the following equation can be 
written down in the absence of nuclear interactions and source contribution: 



d{V{t)N{>E,t)) _ d{V{t)N{>E,t)) dE d{V{t)N{>E,t)) _ 

dt dt dt dE ^ ^ 

that gives 

^?iX^,3//Wm*)-^M^^0, (2) 

where H is the Hubble parameter and the expressions ^ = 3HV and 
^ = —HE are used (the last equation describes the energy loss rate of 
ultrarelativistic particle with energy E ^ pc in the expanding universe), see 
Zeldovich & Novikov (1983). 

Additional terms that describe the continues energy losses -^-miim (if 
an individual particle loses energy as ^ = —7(^7)) and the "catastrophic" 
particle losses of the form v{E,t)N (that describes nuclear fragmentation or 
radioactive decay) can be added to the left side of Eq. (2). The source term 
can be added to its right side. Also, the redshift z can be introduced instead 
of t variable using the relation ^ = —(1 + 2;) i7( 2). Finally, the equation for 
nuclei with mass number A takes the form: 



-H{z){l + z)- 
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(3) 



i=l,2.. 
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The system of equations (3) for all kinds of nuclei with different A should 
be solved simultaneously. The energy per nucleon e = E/A is used here 
because it is approximately conserved in a process of nuclear photodisinte- 
gration (analogously to a standard practice employed in the studies of cosmic 
ray nuclear spallation in the interstellar gas, see Berezinskii et al. (1990)), 
F{A, e, z) is the corresponding cosmic-ray distribution function, q{A, e) is 
the density of cosmic-ray sources at the present epoch 2; = 0, m characterizes 
the source evolution (the evolution is absent for m = 0), T{A^e,z) is the 
characteristic time of energy loss by the production of e~e^ pairs and pions, 
v{A, e, z) is the frequency of nuclear photodisintegration, the sum in the right 
side of Eq. (3) describes the contribution of secondary nuclei produced by 
the photodisintegration of heavier nuclei, H{z) = Hq{{1 + zY^lm + ^a)^''^ is 
the Hubble parameter in the fiat expanding universe with the matter density 
nm{= 0.3) and the A-term fiA(= 0.7). 

We used Eq. (3) for calculations of cosmic-ray flux produced by galaxy 
cluster accretion shocks (Ptuskin et al., 2009) and earlier for calculations 
of proton fluxes generated by different distributions of cosmic ray sources 
(Ptuskin et al., 2003). The analogous equations were recently used by (Aloisio 
et al., 2010). The analytic study of the process of transformation of cosmic- 
ray composition described in (3) by the terms with frequencies u was con- 
ducted by Hooper et al. (2008). They reproduced some mathematical results 
known in the theory of nuclear transformation caused by the spallation reac- 
tions in the interstellar g Berezinskii et al. (1990). Notice that usually 
the authors prefer to build up an expression for energetic particle or photon 
density in the expanding universe without explicit writing down equations 
like (3), see e.g. basic monographs Stecker (1971); Berezinskii et al. (1990). 
Excellent discussion on the description of cosmic ray propagation was pre- 
sented by Berezinsky et al. (2006). 

Eq. (3) are valid for an arbitrary regime of cosmic ray propagation - diffu- 
sion, rectilinear motion, or any intermediate regime. In the case of rectilinear 
propagation, all functions in Eq. (3) may have the direction of particle mo- 
tion (determined by the unit vector — ) as an additional independent variable 
that allows studying some models with anisotropic source distribution. 

The photodisintegration by the background radiation is calculated in the 
present work in the approximation used by Karakula & Tkaczyk (1993). The 
corresponding cross sections are mainly taken from the works Puget et al. 
(1976); Rachen (1996); Khan et al. (2005). The spectra of the background 
microwave, infrared and optical radiation are taken from Malkan & Stecker 
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(1998). The numerical solution of the cosmic-ray transport equations follows 
the finite differences method. The variables are the redshift z and log(£'/yl). 
The computation starts from some maximum z {zmax = 2 was assumed 
in our calculations) and goes in the direction of decreasing redshift till the 
present epoch z = 0. At any given z, the computation starts from some 
maximum energy per nucleon and goes to smaller energies and to smaller 
atomic numbers A from Iron to Hydrogen. 

It should be emphasized that an alternative Monte Carlo techniques were 
used for treating the photodisintegration of ultra-high energy cosmic rays. 
A representative list of references was given by Hooper et al. (2008), see 
also Allard et al. (2008); Allard (2009); Hooper & Taylor (2010) where the 
important results on expected cosmic ray spectra for different source spectra 
and compositions were obtained. Although more sophisticated in treating 
details of nuclear interactions, the Monte Carlo techniques are more time 
consuming compared to the solution of Eq.(3) by the finite differences method 
and should be used when it is required by the specific characteristics of the 
problem under consideration. 

The assumption of continuous source distribution is not valid when parti- 
cles lose energy at a scale less than the distance between cosmic ray sources. 
The finite distance to the nearest source is approximately taken into ac- 
count in our calculations by the cutoff of the source distribution at Zmin ~ 
OASHqcI/c « 1, so that g = at 2; < Zmin {OASd is the average distance of 
an observer to the nearest source if the point sources arrange the cubic lattice 
with the edge, the distance between sources, equals to d). The statistically 
uniform source distribution is assumed at larger redshifts. The rectilinear 
propagation of cosmic rays from the close source is accepted here but dif- 
ferent modes of cosmic ray propagation can be modelled by changing the 
relation between Zmin and d. 

The imposition of an arbitrary intergalactic magnetic field does not affect 
the uniform isotropic distribution of cosmic rays. The account of cosmic-ray 
source granulation makes important the presence of extragalactic magnetic 
fields. Magnetic effects depend on the gyro-radius = 10(-E/10^^ZeV)(i?/10~^G)~^ 
Mpc that determines the deflection of ultrarelativistic particle. In principle, 
the presence of magnetic field changes angular distribution of cosmic rays 
and cause cosmic rays to propagate over longer distances than the straight 
line distance. The weak influence of intergalactic fields on the cosmic ray 
intensity was found in MHD simulations of large-scale structures by Dolag et 
al. (2004). Different conclusion was reached by Sigl et al. (2004). It is clear 
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that results critically depend on the assumptions about the strength and 
structure of the intergalactic magnetic field, see also (Parizot, 2004; Kotera 
& Lemoine, 2008; Das et al., 2008). Some of magnetic field effects can by 
investigated in the diffusion approximation (Berezinsky et al., 2006; Berezin- 
sky & Gazizov, 2007) (additional diffusion terms should be added then to 
the transport equations). The most straightforward but time consuming is 
the direct calculation of possible trajectories in extragalactic magnetic fields 
from a source to an observer. The corresponding Monte Carlo code by Sigl 
et al. (2004) is available on the web. It is worth to emphasize again that the 
main problem in such kind calculations is the uncertainties in the strength 
and structure of extragalactic magnetic fields. 

The objective of our paper, the determination of the basic overall shapes of 
cosmic ray spectrum at the Earth for different source functions produced by 
the AGN jets, does not require more complicated computational procedure 
than solution of Eqs (3) with a source distribution cutoff at some Zmin- 

3. Spectrum of ultra-high energy cosmic rays accelerated in AGN 
jets 

Simple estimates (Torres & Anchordoqui, 2004; Berezinskii et al., 1990) 
show that from the viewpoint of energetics the AGN jets can be the sources 
of ultra-high energy cosmic rays. Different aspects of particle acceleration in 
AGN jets were considered by Biermann & Strittmatter (1987); Norman et 
al. (1995); Biermann et al. (2008); Lemoine (2008); Pe'er et al. (2009), see 
also references below. 

To maintain the cosmic ray intensity observed in the Auger experiment 
at energies above 10^^ eV, the power of extragalactic sources of the order of 
3 X 10^^ erg s~^ Mpc~^ is required. This value increases if the contribution 
of cosmic rays with smaller energies is taken into account. At the same 
time the AGN jets release kinetic energy at the level of 3 x lO^'' erg s~^ 
Mpc~'^ and approximately 6% of this energy is contained in the jets with a 
power Ljet = 10^^ — 10^^ erg s~^ characteristic of FRII (Fanaroff- Riley II) 
radiogalaxies and radio loud quasars. We shall use the notation FII for this 
population of jets. The numerous and less powerful jets of low-luminosty 
AGN have power Lj^t = 10^^ — 10^^ erg s~^. We shall denote them as FI 
sources. 

Without specifying the mechanism of particle acceleration in jets, one 
can use the Hillas criterion Emax = Ze[5Bl for the estimate of maximum 
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energy which the particles with charge Ze in the acceleration region of size 
/, magnetic field strength 5, and the velocity of magnetic field transport 
u = (5c can gain (Hillas (1984), see also Ptitsyna & Troitsky (2010)). 

Let us consider the "optimistic" estimate and assume that the energy fiux 
of a statistically isotropic magnetic field frozen in the jet is related to the 
kinetic energy fiux by the relation Lj^t = Pc^ttR^. The equality R = 1/2 
is accepted here for the jet radius; the jet velocity is /3c. As a result, the 
following estimate of the maximum energy of accelerated particles can be 
obtained: 

E^a. = Ze {QPc-^L^et) ^ 2.7 X lQ''^ZP^I''Lj,tA^^/^eV, (4) 

where Ljet,4,5 = Ljej(10^^erg s~^)~^, see Lovelace (1976); Blandford (1993); 
Aharonian et al. (2002); Waxman (2004); Farrar & Gruzinov (2009) and 
references therein for the derivation of similar formulas. 

The expression for Emax can be also derived based on the well studied 
case of the diffusive shock acceleration in young supernova remnants. Let 
us consider a jet which consists of the proton-electron plasma with the mass 
density p and the power Ljet = O-bpu^nR^ . The cosmic rays are accelerated 
at the jet termination shock and their energy density is Wcr = VcrPu'^, where 
r/cr ~ 0.1. The magnetic field at the site of acceleration can reach the value 
of the order B = (47r/3wcr)^^^ if the field is amplified by the strong cosmic- 
ray streaming instability (Bell, 2004). The maximum energy of accelerated 
particles is E^ax = Ze^BR if the Bohm diffusion near the shock is assumed 
(note that this Emax satisfies the Hillas criterium). It finally gives Emax ~ 
Zef5 {SrjcrC'^Ljett^'^ that is close to the estimate (4). We shall use Eq. (4) 
and set /3 = 1 in the calculations below. 

We consider two types of cosmic-ray source spectrum ejected into the 
intergalactic space. The first type is a delta function spectrum and the cor- 
responding source power is qd = ^crnjetLjetEmax~^S{E - Emax), where the 
coefficient ^cr characterizes the fraction of jet kinetic energy that goes to the 
accelerated particles; rijet is the jet number density in the intergalactic space. 
The second type of sources has a power law spectrum oc E'"^ and the corre- 
sponding source power is = C,crnjetLjetE~'^B.{Emax ~ E), where H(a;) is the 
step function. (Strictly speaking, the additional logarithmic normalization 
factor (\n{Emax/Emin))~^ should by included in the last expression for Qp. We 
omit it because of uncertain value of the minimal energy Emin for particles 
ejected from the accelerator.) 
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Figure 1: Calculated average source spectra of jet populations FI (solid lines) and FII 
(dash line) for delta- function (d) and power law (p) cosmic ray spectra generated by 
individual jets. Eq. (4) was used for calculations of Emax{Ljet) 

It can be recalled as an example that both spectrum shapes of ejected 
particles arise in the consideration of diffusive shock acceleration in supernova 
remnants, see Ptuskin & Zirakashvili (2005). The runaway particles have 
close to the delta function energy spectrum ~ 5{E — Emax) where Emax is 
the maximum energy of accelerated particles that is achieved at the given 
stage of a supernova shock evolution. The energetic particles that remains 
confined inside the remnant may have close to a power law spectrum and 
leave out into the interstellar medium at some stage of SNR evolution when 
the shock breaks up. It should be stressed that we do not assume that 
two discussed injection spectra work at the same time and analyze them 
separately. 

The source functions qd{E) and qp{E) should be averaged over the distri- 
bution of jet luminosity njet{Ljet) to obtain the average source function of 
extragalactic cosmic rays. The results are shown in Fig. (1) for the func- 
tions njet{Ljet) presented by Koerding et al. (2008) in their Figure (8) for 
the kinetic luminosity function of jets. Four source functions in Fig. (1) cor- 
respond to the combination of two populations of jets, FI and FII, and two 
types of jet spectra, d and p. These four source functions are considered here 
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Figure 2: Calculated spectra of extragalactic cosmic rays for sources < qd > and < qp > 
averaged over the AGN jet population FI (solid lines Fid and FIp) and FII (dash lines 
Flld and Flip). The spectra (except Flld) are normalized to the intensity observed at 
10^^ eV in the Auger experiment. Data are from HiRes experiment (Abbasi et al., 2009) 
(circles) and Auger experiment (Abraham et al., 2010) (squares). 

as representing four different scenarios of cosmic ray acceleration in jets. 

We solved numerically the set of transport Eqs (3) for the described types 
of the average source spectrum. Fig. (2) illustrates the results. It was as- 
sumed that the source chemical composition coincides with the composition 
of Galactic cosmic ray sources. The spectra were normalized at 10^^ eV to 
the observed by Auger intensity. It requires very different efficiency of par- 
ticle acceleration ^cr,Fii/^cr,Fi ~ 20 in the FII and FI jets. These efficiencies 
are ^cr,Fii ~ 0.1 and xicr,Fi ~ 0.005 if the source spectra are extrapolated 
down to 1 GeV. No cosmological evolution was assumed in our calculations 
(m = 0). The evolution is different for different morphological types of AGN, 
see e.g. Berezinsky et al. (2006), but it does not significantly affect the calcu- 
lated spectra at energies > 3 x 10^^ eV since these particles may come from 
the distances not larger than about 2 x 10^ Mpc. 

Of four spectra shown in Fig. (2), two reproduce cosmic ray observations 
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Figure 3: Calculated spectrum of extragalactic cosmic rays for sources < qd > averaged 
over the AGN jet population Flld with heavy composition; Emax is decreased by 80 
compared to Eq. (4). Data are from Auger experiment (Abraham et al., 2010) (squares). 



with reasonable accuracy. They correspond to the AGN jet population FI 
with delta-function source spectra and the AGN jet population FII with 
power-law spectra E'"^. 

The finite distance Zmin to the closest to an observer source was taken into 
account as discussed in the preceding Section. This distance is a function of 
particle energy and charge and is different for the source populations FI and 
FII. The absence of sources at distances < 90 Mpc in our model resulted in 
a steeply sloping down spectrum of cosmic rays at the highest energies 10^° 
eV for the FII source distribution. 

The dependence njet{Ljet) together with Eq. (4) for Emax{Ljet) leads to the 
dependence of cosmic-ray source number density on particle energy ns{E). 
For FI population of sources with delta shaped spectra, the source density is 
Us = 10"^ Mpc-3 at E = 6x 10^^ eV and = 2 x lO'^ Mpc'^ at E = 10^^ 
eV that coincides with the results of Takami & Sato (2009) derived from the 
analysis of cosmic-ray arrival direction distribution. 

Protons dominate in the calculated composition of cosmic rays that is in 
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strong disagreement with the Auger data. To get out of a difficulty one can 
take anomalously high abundance of heavy nuclei and reduce the maximum 
particle energy at the source (Allard, 2009; Hooper & Taylor, 2010). We re- 
serve the consideration of this issue for a later paper and show in Fig. (3) only 
one example of calculations where the shape of source spectrum corresponds 
to the Flld sources with the maximum particle energy Emax decreased by 
a factor of 80 compared to the "optimistic" estimate (4) (this relieves the 
extreme assumptions used in the derivations of Eq. (4)) and the Iron abun- 
dance at the source comprises 1/3 of all nuclei. The calculated cosmic-ray 
spectrum only roughly reproduce the observed spectrum. The cosmic-ray 
composition can be characterized by the mean logarithmic atomic number 
< ln(A) >. Its calculated value rises approximately linearly from about 0.25 
at 5 X 10^* eV to 3.7 at 5 x 10^^ eV in a qualitative agreement with the Auger 
data on the shower maximum dependence on energy. 

4. Conclusions 

It is believed that the particle acceleration by jets in active galactic nuclei 
is the most efficient source of cosmic rays with the highest energies, E > 10^^ 
eV. In the present work we distinguish two populations of jets: FI produced 
by the low luminosity AGN populations with jet power 2 x 10"^° to 3 x 10^^ 
erg s^^, and FII produced by high luminosity AGN with jet power larger 
than 2 x 10^^. The corresponding jet distributions on power were given by 
Koerding et al. (2008). 

The typical power law spectrum of nonthermal jet radiation implies the 
power law particle spectrum of the form close to E'"^. One may expect that 
the spectrum of particles released into the extragalactic space, the source 
spectrum of extragalactic cosmic rays, has the same shape. Another pos- 
sibility is that accelerated particles remains confined inside the source and 
only particles with maximum energies run away into intergalactic space so 
that the source spectrum is of a delta-shaped form. The average source spec- 
trum of extragalactic cosmic rays is determined as the convolution of one of 
these source functions of an individual jet with the jet distribution on power. 
Based on the Hillas criterion, we accepted an optimistic estimate for the 
maximum energy of accelerated particles (4) with its characteristic scaling 
Emax oc and used it in our calculations. 

The computations of cosmic ray propagation in the expanding Universe 
filled with the background electromagnetic radiation were fulffiled by the 
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use of a simple numerical code which solves the system of coupled transport 
Eqs (3) for energetic protons and nuclei from He to Fe. The calculations were 
made under the approximation of continuous energy losses by e~, e+ and 
pion production and the "catastrophic" losses through photodisintegration 
and corresponding production of secondary nuclei. 

The results of our calculations are illustrated in Fig. (2). The observed 
spectrum of ultra-high energy cosmic rays can in principle be explained in 
the frameworks of two scenarios - the acceleration by the FI sources with 
individual jet spectra of the delta-shaped form or the acceleration by the FIT 
sources with individual jet spectra close to E""^ form. The transition from 
Galactic to extragalactic component in the observed at the Earth spectrum 
occurs at about (3. ..5) x lO^^eV. 

The calculated spectra are normalized to the observed by Auger inten- 
sity at 10^^ eV. It requires very different efficiencies of transformation of 
the jet kinetic energy to the energy of cosmic rays in the FII and FI sources: 
^cr,Fii/^cr,Fi ~ 20. It was assumcd that the elemental composition of acceler- 
ated particles is the same as in the Galactic cosmic-ray sources. This results 
in the predominantly proton composition of ultra high energy extragalactic 
cosmic rays that is compatible with the HiRes data but not with the Auger 
data. One needs to significantly increase the abundance of heavy nuclei at 
the source and drastically decrease the value of E^ax to fit the Auger data 
(Allard, 2009; Hooper & Taylor, 2010). This procedure was discussed at the 
end of Section 3 and illustrated by Fig. 3. 

The work was supported by the Russian Foundation for Basic Research 
grant 10-02-00 110a. 
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